
require(Rcpp)
options(rasterTmpDir='C:/temp/')
output_version <- '2020-09-14' # Sys.Date()
  
ntl.zones <- read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_ntl_adm0_x_d01') %>%
  st_set_crs(.,4326)
names(ntl.zones)
if ((length(unique(ntl.zones$OBJECTID))==dim(ntl.zones)[1])==FALSE) { break }

##
## Built-up GHS
##
## https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/GHSL/GHS_BUILT_LDSMT_GLOBE_R2018A
ghs_fn = list.files(ghs_dir, pattern=".tif$", recursive = TRUE, full.names=TRUE)
ghs_stack = stack(ghs_fn)
built.zones.moll <- st_transform(subset(ntl.zones,select=c(OBJECTID)),crs(ghs_stack))
built.zones.moll$area_m2 <- st_area(built.zones.moll)
# crop to study area
ghs_built_crop = crop(ghs_stack,built.zones.moll,snap='out')

# Built-up area density by epoch, aggregated from 30 m.
# 250m of resolution - World Mollweide (EPSG:54009)
# share of built up * res
ghs_built_m2 = ghs_built_crop * (res(ghs_stack)[1]) * (res(ghs_stack)[2])
built.extr <- exact_extract(ghs_built_m2,built.zones.moll,fun='sum',include_xy=FALSE, weights = NULL, progress=TRUE,max_cells_in_memory = 3e+07)

##
## GHS-POP
##

##
ghspop_fn = list.files(ghspop_dir, pattern=".tif$", recursive = TRUE, full.names=TRUE)
ghspop_stack = stack(ghspop_fn)
pop.zones.moll <- st_transform(subset(ntl.zones,select=OBJECTID),crs(ghspop_stack))
# crop to study area
ghs_pop_crop = crop(ghspop_stack,pop.zones.moll,snap='out')
pop.extr <- exact_extract(ghs_pop_crop,pop.zones.moll,fun='sum',include_xy=FALSE, weights = NULL, progress=TRUE,max_cells_in_memory = 3e+07)
ntl.dt = dplyr::bind_cols(ntl.zones,built.extr,st_drop_geometry(built.zones.moll)["area_m2"],pop.extr)
names(ntl.dt)

df <- as.data.frame(ntl.dt)
#names(df[,c(2,4,33:length(names(ntl.dt)))])
df1 <- df %>%
  dplyr::select(OBJECTID,GID_0,
                sum.GHS_BUILT_LDS1975_GLOBE_R2018A_54009_250_V2_0,
                sum.GHS_BUILT_LDS1990_GLOBE_R2018A_54009_250_V2_0,
                sum.GHS_BUILT_LDS2000_GLOBE_R2018A_54009_250_V2_0,
                sum.GHS_BUILT_LDS2014_GLOBE_R2018A_54009_250_V2_0,
                area_m2,
                sum.GHS_POP_E1975_GLOBE_R2019A_54009_250_V1_0,
                sum.GHS_POP_E1990_GLOBE_R2019A_54009_250_V1_0,
                sum.GHS_POP_E2000_GLOBE_R2019A_54009_250_V1_0,
                sum.GHS_POP_E2015_GLOBE_R2019A_54009_250_V1_0
                )

names(df1) <- c("objectid","gid_0",paste0("built",c(1975,1990,2000,2015)),"area_m2",paste0("pop",c(1975,1990,2000,2015)))
save.dta13(df1, paste0(wkdir,"/proc_data/ghs_built_pop",output_version,".dta"))